function [beta,lambda0,lambda1,tstatbeta,selambda0,selambda1,r2]=adrianetal(f,r,nc,nf)
% r return vector T*N
% f state factor vector T*K
% kc columns of f that are priced
% kf columns of f that drive time variation of risk prices
[T,~] = size(f); % sample length of the state factors
[~,N] = size(r);
% Step 1 VAR estimate
X = [ones(T-1,1) f(1:end-1,:)];
Y = f(2:end,:);
Phi =(X'*X)\(X'*Y);
u = Y-X*Phi; % residuals from VAR
u_c = u(:,nc);
sigmau = u_c'*u_c/(T-1);
X_f = [X(:,1) f(1:end-1,nf)];
sigmaf = X_f'*X_f/(T-1);
% Step 2 Reduced-form SUR regression
u_f = u(:,nc); % Extract the columns of risk factors
t_f = f(1:end-1,nf); % Extract the columns of price of risk factors
r = r(2:end,:);
X2 = [t_f u_f];
X3 = [ones(T-1,1) t_f u_f];

kc = length(nc);
kf = length(nf);

a0 = nan(N,1);
a1 = nan(kf,N);
beta = nan(kc,N);
tstatbeta = nan(kc,N);
r2 = nan(N,1);
e = nan(T-1,N);

ind_r = isnan(r);
sum_ind_r = sum(ind_r,2);
start = min(find(sum_ind_r==0));

for i = 1:N
result = fitlm(X2,r(:,i));
a0(i,1) = result.Coefficients.Estimate(1,1);
a1(:,i) = result.Coefficients.Estimate(2:kf+1,1);
beta(:,i) = result.Coefficients.Estimate(kf+2:end,1);
tstatbeta(:,i) = result.Coefficients.tStat(kf+2:end,1);
r2(i,1) = result.Rsquared.Ordinary;
e(:,i) = result.Residuals.Raw(:,1);
end
a1 = a1';
beta = beta'; %N*kc
tstatbeta = tstatbeta';
% Heteroskedasticity robust variance variance matrix
X3_start = X3(start:end,:);
T_start = length(X3_start(:,1));
e_start = e(start:end,:);

vrob_1 = kron(inv(X3_start'*X3_start),eye(N));
vrob_2 = zeros(N*(kc+kf+1),N*(kc+kf+1));
for i = 1:T_start
    temp = kron(X3_start(i,:)'*X3_start(i,:),e_start(i,:)'*e_start(i,:));
    vrob_2 = vrob_2+temp;
end
vrob = T_start*vrob_1*vrob_2*vrob_1;

% Step 3
lambda0 = (beta'*beta)\(beta'*a0);
if kf>0
lambda1 = (beta'*beta)\(beta'*a1);
end
if kf == 0
    lambda1 = [];
end
lambda = [lambda0 lambda1];
% Standard error estimate
h1 = kron(eye(kf+1),(beta'*beta)\beta');
h2 = -kron(lambda',(beta'*beta)\beta');
h = [h1 h2];
vlambda = 1/T_start*(kron(inv(sigmaf),sigmau)+h*vrob*h');
selambda0 = sqrt(diag(vlambda(1:kc,1:kc)));
selambda1 = sqrt(diag(vlambda(kc+1:end,kc+1:end)));
if kf>0
lambda1 = reshape(lambda1,kc*kf,1);
end
end